#include <stdio.h>
#include <math.h>
#include "bmath.h"

double gammafn(double x) {
    const static double gamcs[42] = {
        +.8571195590989331421920062399942e-2,
        +.4415381324841006757191315771652e-2,
        +.5685043681599363378632664588789e-1,
        -.4219835396418560501012500186624e-2,
        +.1326808181212460220584006796352e-2,
        -.1893024529798880432523947023886e-3,
        +.3606925327441245256578082217225e-4,
        -.6056761904460864218485548290365e-5,
        +.1055829546302283344731823509093e-5,
        -.1811967365542384048291855891166e-6,
        +.3117724964715322277790254593169e-7,
        -.5354219639019687140874081024347e-8,
        +.9193275519859588946887786825940e-9,
        -.1577941280288339761767423273953e-9,
        +.2707980622934954543266540433089e-10,
        -.4646818653825730144081661058933e-11,
        +.7973350192007419656460767175359e-12,
        -.1368078209830916025799499172309e-12,
        +.2347319486563800657233471771688e-13,
        -.4027432614949066932766570534699e-14,
        +.6910051747372100912138336975257e-15,
        -.1185584500221992907052387126192e-15,
        +.2034148542496373955201026051932e-16,
        -.3490054341717405849274012949108e-17,
        +.5987993856485305567135051066026e-18,
        -.1027378057872228074490069778431e-18,
        +.1762702816060529824942759660748e-19,
        -.3024320653735306260958772112042e-20,
        +.5188914660218397839717833550506e-21,
        -.8902770842456576692449251601066e-22,
        +.1527474068493342602274596891306e-22,
        -.2620731256187362900257328332799e-23,
        +.4496464047830538670331046570666e-24,
        -.7714712731336877911703901525333e-25,
        +.1323635453126044036486572714666e-25,
        -.2270999412942928816702313813333e-26,
        +.3896418998003991449320816639999e-27,
        -.6685198115125953327792127999999e-28,
        +.1146998663140024384347613866666e-28,
        -.1967938586345134677295103999999e-29,
        +.3376448816585338090334890666666e-30,
        -.5793070335782135784625493333333e-31
    };

    int i;
    double y;
    double value;

#ifdef NOMORE_FOR_THREADS
    static int ngam = 0;
    static double xmin = 0, xmax = 0., xsml = 0., dxrel = 0.;

    /* Initialize machine dependent constants, the first time gamma() is called.
        FIXME for threads ! */
    if (ngam == 0) {
        ngam = chebyshev_init(gamcs, 42, DBL_EPSILON / 20); /*was .1*d1mach(3)*/
        gammalims(&xmin, &xmax); /*-> ./gammalims.c */
        xsml = exp(fmax2(log(DBL_MIN), -log(DBL_MAX)) + 0.01);
        /*   = exp(.01)*DBL_MIN = 2.247e-308 for IEEE */
        dxrel = sqrt(DBL_EPSILON); /*was sqrt(d1mach(4)) */
    }
#else
    /* For IEEE double precision DBL_EPSILON = 2^-52 = 2.220446049250313e-16 :
     * (xmin, xmax) are non-trivial, see ./gammalims.c
     * xsml = exp(.01)*DBL_MIN
     * dxrel = sqrt(DBL_EPSILON) = 2 ^ -26
     */
#define ngam 22
#define xmin -170.5674972726612
#define xmax  171.61447887182298
#define xsml 2.2474362225598545e-308
#define dxrel 1.490116119384765696e-8
#endif

    if (isnan(x)) return x;

    /* If the argument is exactly zero or a negative integer
     * then return NaN. */
    if (x == 0 || (x < 0 && x == (long) x)) {
        return NAN;
    }

    y = fabs(x);

    if (y <= 10) {

        /* Compute gamma(x) for -10 <= x <= 10
         * Reduce the interval and find gamma(1 + y) for 0 <= y < 1
         * first of all. */

        int n = floor(x);
        if (x < 0) --n;
        y = x - n; /* n = floor(x)  ==>	y in [ 0, 1 ) */
        --n;
        value = chebyshev_eval(y * 2 - 1, gamcs, ngam) + .9375;
        if (n == 0)
            return value; /* x = 1.dddd = 1+y */

        if (n < 0) {
            /* compute gamma(x) for -10 <= x < 1 */

            /* exact 0 or "-n" checked already above */

            /* The argument is so close to 0 that the result would overflow. */
            if (y < xsml) {
                if (x > 0) return INFINITY;
                else return -INFINITY;
            }

            n = -n;

            for (i = 0; i < n; i++) {
                value /= (fabs(x + i) > 10.0e-15) ? (x + i): 1.0;
            }
            return value;
        } else {
            /* gamma(x) for 2 <= x <= 10 */

            for (i = 1; i <= n; i++) {
                value *= (y + i);
            }
            return value;
        }
    } else {
        /* gamma(x) for	 y = |x| > 10. */

        if (x > xmax) { /* Overflow */
            return INFINITY;
        }

        if (x < xmin) { /* Underflow */
            return 0.;
        }

        if (y <= 50 && y == (int) y) { /* compute (n - 1)! */
            value = 1.;
            for (i = 2; i < y; i++) value *= i;
        } else { /* normal case */
            value = exp((y - 0.5) * log(y) - y + M_LN_SQRT_2PI +
                    ((2 * y == (int) 2 * y) ? stirlerr(y) : lgammacor(y)));
        }
        if (x > 0)
            return value;

        double sinpiy = sin(M_PI * y);
        if (sinpiy == 0) { /* Negative integer arg - overflow */
            return INFINITY;
        }

        return -M_PI / (y * sinpiy * value);
    }
}